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Abstract 

Classical effective potentials are indispensable for any large-scale atomistic sim- 
ulations, and the relevance of simulation results crucially depends on the quality of 
the potentials used. For complex alloys like quasicrystals, however, realistic effective 
potentials are practically inexistent. We report here on our efforts to develop effec- 
tive potentials especially for quasicrystalline alloy systems. We use the so-called force 
matching method, in which the potential parameters are adapted so as to optimally 
reproduce the forces and energies in a set of suitably chosen reference configurations. 
These reference data are calculated with ab-initio methods. As a first application, 
EAM potentials for decagonal Al-Ni-Co, icosahedral Ca-Cd, and both icosahedral 
and decagonal Mg-Zn quasicrystals have been constructed. The influence of the po- 
tential range and degree of specialisation on the accuracy and other properties is 
discussed and compared. 

Keywords: force matching; quasicrystal; effective potential; molecular dynamics; ab 
initio 



1 Introduction 

Large-scale molecular dynamics simulations are possible only with classical effective po- 
tentials, which reduce the quantum-mechanical interactions of electrons and nuclei to an 
effective interaction between the atom cores. The computational task is thereby greatly 
simplified. Whereas ab-initio simulations are limited to a few hundred atoms at most, 
classical simulations can be done routinely with multi-million atom systems. For many 
purposes, such system sizes are indispensable. For example, fracture studies of quasicrys- 
tals require samples with several million atoms at least [1]. Diffusion studies, on the other 
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hand, can be done with a few thousand atoms (or even less), but require very large simu- 
lated times of the order of nanoseconds [2], which also makes them infeasible for ab-initio 
simulations. 

While physically justified effective potentials have been constructed for many elemen- 
tary solids, such potentials are rare for complex intermetallic alloys. For this reason, molec- 
ular dynamics simulations of these materials have often been done with simple model po- 
tentials, resulting in rather limited reliability and predictability. In order to make progress, 
better potentials are needed to accurately simulate complex materials. 

The force matching method [3] provides a way to construct physically reasonable po- 
tentials also for more complex solids, where a larger variety of local environments has to 
be described correctly, and many potential parameters need to be determined. The idea is 
to compute forces and energies from first principles for a suitable selection of small refer- 
ence systems, and to fit the potential parameters so that they optimally reproduce these 
reference data. Hereafter, potentials generated in this way will be referred to as fitted 
potentials. Thus, the force matching method allows to make use of the results of ab-initio 
simulations also for large-scale classical simulations, thereby bridging the gap between the 
sample sizes supported by these two methods. 

2 Force Matching 

As we intend to construct potentials for complex intermetallic alloys, we have to assume 
a functional form which is suitable for metals. A good choice are EAM (Embedded Atom 
Method) potentials [4], also known as glue potentials [5]. Such potentials have been used 
very successfully for many metals, and are still efficient to compute, even though they 
include many-body terms. In contrast, pure pair potentials show a number of deficiencies 
when it comes to describe metals [5] . The functional form of EAM potentials is given by 

V = ^Mi( r «) + U ^( U i)> with U i = ^Pkjinj), (i) 

i,j<i i j^i 

where 0*^. is a pair potential term depending on the two atom types k\. describes 
the embedding term that represents an additional energy for each atom. This energy is a 
function of a local density n« determined by contributions pk of the neighbouring atoms. 
It is tempting to view this as embedding each atom into the electron sea provided by 
its neighbours. Such an interpretation is not really meaningful, however. The potential 
([!]) is invariant under a family of "gauge" transformations [5], by which one can move 
contributions from the embedding term to the pair term, and vice versa, so that it makes 
little sense to give any of them an individual physical interpretation. 

In order to allow for maximal flexibility, and to avoid any bias, the potential functions 
in ([!]) are represented by tabulated values and spline interpolation, the tabulated values 
acting as potential parameters. This makes it unnecessary to guess the right analytic form 
beforehand. The sampling points can be chosen freely, which is useful for functions which 
vary rapidly in one region, but only slowly in another region. 
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The forces and energies in the reference structures are computed with VASP, the Vi- 
enna Ab-initio Simulation Package [6,7], using the Projector Augmented Wave (PAW) 
method [8,9]. Like all plane wave based ab-initio codes, VASP requires periodic bound- 
ary conditions. For quasicrystals, this means that periodic approximants have to be used 
as reference structures. As ab-initio methods are limited to a few hundred atoms, those 
approximants must be rather small. For the systems studied so far, this was not a major 
problem, as the relevant local environments in the quasicrystal all occur also in reasonably 
small approximants. Icosahedral quasicrystals with F-type lattice may be more problem- 
atic in this respect. For these, small approximants are rare, and the force matching method 
requires a sufficient variety of reference structures. 

Given the reference data (forces, energies, and stresses in the reference structures), 
the potential parameters (in our case: up to about 120 EAM potential sampling points 
for spline interpolation) then are optimised in a non-linear least square fit, so that the 
fitted potential reproduces the reference data as well as possible. The target function to be 
minimised is a weighted sum of the squared deviations between the reference data, denoted 
by the subscript below, and the corresponding data computed from the fitted effective 
potential. It is of the form 

Z = Z F + Z C , with (2) 
* = E E W ^;- f ^ )2 , and Z c = f:W k ^-y\ (3) 

j=l a=x,y,z J 0,j + E 3 k=1 ^O.fc + E k 

where Zp represents the contributions of the forces / -, and Zc those of some collective 
quantities like total stresses and energies, but also additional constraints A k on the potential 
one would like to impose. The denominators of the fractions ensure that the target function 
measures the relative deviations from the reference data, except for really tiny quantities, 
where the ei prevent extremely small denominators. The W\ are the weights of the different 
terms. It proves useful for the fitting to give the total stresses and the cohesion energies a 
higher weight, although in principle they should be reproduced correctly already from the 
forces. 

We developed a programme named potfit, which optimises the potential parameters to 
a set of reference data. It consists of two largely independent parts. The first part imple- 
ments a particular parametrised potential model. It takes a list of potential parameters 
and computes from it the target function, i.e., the deviations of the forces, energies, and 
stresses from the reference data. Wrapped around this part is a second, potential indepen- 
dent part, which implements a least square minimisation module, using a combination of a 
deterministic conjugate gradient algorithm [10] and a stochastic simulated annealing algo- 
rithm [11]. This part knows nothing about the details of the potential, and only deals with 
a list of potential parameters. The programme architecture thus makes it easy to replace 
the potential dependent part by a different one, e.g., one which implements a different 
potential model, or a different way to parametrise it. 
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3 Results and Applications 



We generated several fitted potentials for decagonal Al-Ni-Co and icosahedral Ca-Cd qua- 
sicrystals, as well as Mg-Zn potentials suitable for both icosahedral and decagonal phases. 
In a first step, classical molecular dynamics simulations with simple model potentials were 
used to create reference configurations from small approximants (80-250 atoms). These in- 
cluded samples at different temperatures, but also samples which were scaled and strained 
in different ways. The approximants were carefully selected, so that all relevant local envi- 
ronments are represented. For those reference structures, the forces, stresses and energies 
were computed with ab-initio methods, and a first version of the fitted effective potential 
given by sampling points with cubic spline interpolation was fitted to the reference data. 
In a second step, molecular dynamics simulations with the newly determined potential 
were used to create new reference structures, which are better representatives of the struc- 
tures actually appearing in that system. The new reference structures complemented and 
partially replaced the previous ones, and the fitting procedure was repeated. This second 
iteration resulted in a significantly better fit to the reference data. In order to test the 
transferability of the fitted potentials, further samples similar to the reference structures 
were created, and their ab-initio forces and energies were compared to those determined 
by the classical potentials. The deviations were of the same order as the deviations found 
in the potential fit, which shows that the fitted potentials transfer well to similar struc- 
tures. For Al-Ni-Co, a force-matched potential is displayed in figure [T] Fitted potentials 
for Ca-Cd and Mg-Zn are not displayed here for space constraints, but are available from 
the authors. 

The potentials developed for decagonal Al-Ni-Co quasicrystals are intended to be used 
in high-temperature diffusion simulations [2]. It is therefore important that they describe 
high temperature states well, which is achieved by selecting the reference structures ac- 
cordingly. By using high temperature reference structures, the fitted potential is especially 
trained to such situations. As part of the potential validation, the melting temperature was 
determined by slowly heating the sample at constant pressure, and the elastic constants of 
decagonal Al-Ni-Co were determined. We actually have constructed two potential variants: 
Variant A gives excellent values for the elastic constants (Table [T]), but produces a melting 
temperature which is somewhat too high. Conversely, variant B shows larger deviations 
in the elastic constants, but gives a very reasonable value of the melting temperature of 
about 1300 K. It is a general experience that with an effective potential it is often not 
possible to reproduce all desired quantities equally well at the same time. 

In complex intermetallic systems there are many competing candidates for the ground 
state structure. This is the case also for complex crystalline systems. In principle, the 
ground state of these can be determined directly by ab-initio simulations, but for large 
unit cells this is extremely time-consuming, or even impossible. Classical potentials can be 
used to select the most promising candidates, and to pre-relax them, so that the time for 
ab-initio relaxation can be dramatically reduced. Potentials used for this purpose must be 
able to discriminate energy differences of the order of a meV/atom. This has been largely 
achieved with fitted potentials for the Mg-Zn and Ca-Cd systems, by using mainly near 
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Figure 1: Potential functions for decagonal Al-Ni-Co 



ground state structures as reference structures. Also, for this application it is important to 
choose a small Sj in equation (|3]), so that small forces are also reproduced accurately. The so 
constructed Ca-Cd potentials have been used successfully for structure optimisations [13]. 



4 Discussion and Conclusion 

The selection of the reference structures used for the potential fit largely determines the 
capabilities of the resulting potential. For a precise determination of the ground state, 
low temperature structures should be dominant in the reference structures, and it must be 
assured that even small forces and energy differences are reproduced accurately. For high 
temperature simulations, on the other hand, typical high temperature structures must be 
predominant in the reference structures. This opens the possibility to design specialised 
potentials for certain purposes by a suitable selection of reference structures. It should 
be kept in mind, however, that a fitted potential can only deal with situations it has 
been trained to. For instance, one should not expect a fitted potential to handle surfaces 
correctly, if it was trained only with bulk systems. Clearly, there is always a trade-off 
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Table 1: Elastic constants of decagonal Al-Ni-Co 
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between the transferability and the accuracy of a fitted potential. A potential can be 
made more versatile by training it with many different kinds of structures, but the more 
versatile it becomes, the less accurate it will be on average. Conversely, very accurate fitted 
potentials will probably have limited transferability. 

For practical applications, the range of a potential is also an important issue, as it 
enters in the third power in the computational effort of molecular dynamics. Allowing for 
a larger potential range results in greater flexibility of the potential, which might improve 
its accuracy, but this comes at the price of a slower simulation. We therefore need a 
compromise between speed and accuracy. The potential range should only be increased as 
long as this can improve the potential quality. In a first step, our fitted potentials were 
constructed with a fairly generous range of about 7A. It turned out, however, that especially 
the transfer function did not make effective use of this range, and was essentially zero 
beyond 5A. In a second fit we therefore restricted the range of to 5A, without significant 
loss of accuracy. This is one of the advantages of using tabulated functions: The system 
itself chooses the optimal functions, including the optimal range. 

Force Matching has proven to be a versatile method to construct physically reasonable, 
accurate effective potentials even for structures as complicated as quasicrystals and their 
approximants. Our potfit programme makes it easy to apply this method to different sys- 
tems, and is also easy to adapt for the support of further potential models. The potentials 
constructed so far have successfully been used in high temperature diffusion simulations 
of decagonal Al-Ni-Co [2], and in structure optimisation of approximants in the Zn-Mg 
and Ca-Cd systems. Further fruitful applications of the fitted potentials can certainly be 
expected, and we hope to apply our methods also to other complex alloy systems, where 
reliable potentials are still lacking. 
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